function [solution] = reconstruct_branch_fault(G, fault_eid, init_branch_status, BusTable, BranchTable, SubstationTable)
%RECONSTRUCT_BRANCH_FAULT 此处显示有关此函数的摘要
%   支路故障的故障隔离分为4种类别。
%   对于双端均有开关的线路，故障发生时直接切除线路即可,类型标记为1
%   对于开关装设在首端的线路，故障发生时等同于末端节点故障，类型标记为2
%   对于开关装设在末端的线路，故障发生时等同于首端节点故障，类型标记为3
%   对于没有开关的线路，故障发生时等同于任选首末节点故障，类型标记为0

    yalmip('clear')
    solution = [];
    Nbus = numnodes(G);
    Nbr = numedges(G);
    Nsrc = height(SubstationTable);
    swidx = find(G.Edges.fixed == false);
    Nsw = numel(swidx);
    Nfixed = Nbr - Nsw;
    br_ind = [G.Edges(:, 'name'), table((1:Nbr).', 'VariableNames', {'index'})];
    br_ind = innerjoin(br_ind, BranchTable(:, {'name', 'MVA', 'r', 'x','is_sw','sw_type'}), 'Keys', 'name', 'RightVariables', {'MVA', 'r', 'x','is_sw','sw_type'});
    br_ind = sortrows(br_ind, 'index', 'ascend');
    bus_ind = sortrows(BusTable, 'busid', 'ascend');


    %故障隔离
    if br_ind.sw_type(fault_eid)==1
        branch_isolated_table=G.Edges(fault_eid, :);
        branch_isolated_table.eid=fault_eid;
    else
        switch br_ind.sw_type(fault_eid)
            case 2
                [~,fault_vid]=findedge(G,fault_eid);
            case 3
                [fault_vid,~]=findedge(G,fault_eid);
            case 0
                [~,fault_vid]=findedge(G,fault_eid);
        end
        if ~isempty(intersect(fault_vid, SubstationTable.busid))
            error('支路%d故障会导致电源被切除，目前程序没有处理这种状况', fault_eid)
        end

        idx = find(~G.Edges.fixed);
        H = rmedge(G, idx);
        bus_isolated = bfsearch(H, fault_vid, 'Restart', false);
        branch_isolated = [];
    
        for k = 1:length(bus_isolated)
            branch_isolated = union(branch_isolated, outedges(G, bus_isolated(k)));
        end
    
        branch_isolated_table = G.Edges(branch_isolated, :);
        branch_isolated_table.eid = reshape(branch_isolated, [], 1);
    end


    %隔离后供电恢复，最小化失电负荷
    sw_status = binvar(Nsw, 1);
    sw_status_table = table(reshape(swidx, [], 1), reshape(1:Nsw, [], 1), ...
        'VariableNames', {'eid', 'varid'}); %记录开关eid和优化变量的对应关系

    Pij = sdpvar(Nbr, 2);
    is_restored = binvar(Nbus, 1);
    Pi = sdpvar(Nbus, 1);
    Psub = sdpvar(Nsrc, 1);
    %System operation constraints
    operation_constraints = [Pij >= -1.0 * [br_ind.MVA, br_ind.MVA], Pij <= [br_ind.MVA, br_ind.MVA]];
    operation_constraints = [operation_constraints, Psub >= 0, Psub <= SubstationTable.MVA];

    %topology constraints
    %1.隔离故障区域的开关固定为0；2.松弛的辐射状约束
    opensw_table = innerjoin(branch_isolated_table(branch_isolated_table.fixed == false, :), sw_status_table, ...
        'Keys', {'eid'}, 'LeftVariables', {'eid'}, 'RightVariables', {'varid'});
    topology_constraints = [sw_status(opensw_table.varid) == 0, ...
                                sum(sw_status) + Nfixed <= (Nbus - (Nsrc + 1))];

    %system power flow constraints
    powerflow_constraints = (is_restored(SubstationTable.busid) == 1);

    % 网络中的支路分为3类，故障区域支路，非故障区域不可动支路，非故障区域可动支路
    % 1.故障区域支路，支路功率固定为0，
    % 2.非故障区域不可动支路，对支路功率作限制
    % 3.非故障区域可动支路，开关状态为1时，对支路功率作限制，开关状态为0时，支路功率固定为0.
    branch_constraint_type = int32(zeros(Nbr, 1));
    branch_constraint_type(G.Edges.fixed == true) = int32(2);
    branch_constraint_type(G.Edges.fixed == false) = int32(3);
    branch_constraint_type(branch_isolated_table.eid) = int32(1);

    for k = 1:Nbr
        cons = [];

        switch branch_constraint_type(k)
            case 1
                cons = (Pij(k, [1 2]) == 0);
            case 2
                cons = (Pij(k, 1) + Pij(k, 2) == 0);
            case 3
                varid = sw_status_table.varid(sw_status_table.eid == k);
                cons = [ ...
                            implies(sw_status(varid) == 1, Pij(k, 1) + Pij(k, 2) == 0), ...
                            implies(sw_status(varid) == 0, Pij(k, [1 2]) == 0) ...
                        ];
        end

        powerflow_constraints = [powerflow_constraints, cons];
    end

    % 节点功率约束,对所有节点有效
    for k = 1:Nbus
        outBranch = outedges(G, k);
        [s_out, t_out] = findedge(G, outBranch);
        s_idx = ismember(s_out, k);
        s_out = outBranch(s_idx);
        t_idx = ismember(t_out, k);
        t_out = outBranch(t_idx);

        Pi(k) = -1.0 * bus_ind.Pd(k);
        idx = (SubstationTable.busid == k);

        if any(idx)
            Pi(k) = Pi(k) + Psub(idx);
        end

        cons = [ ...
                    implies(is_restored(k) == 1, sum(Pij(s_out, 1)) + sum(Pij(t_out, 2)) == Pi(k)), ...
                    implies(is_restored(k) == 0, [Pij(s_out, 1) == 0, Pij(t_out, 2) == 0]) ...
                ];
        powerflow_constraints = [powerflow_constraints, cons];
    end

    %Object
    %Objective = sum(abs(branch_status - reshape(init_sw_status, [], 1)));%最小化开关次数
    %Objective = -1.0 * sum(is_restored .* bus_ind.Pd); %最大化恢复负荷
    init_sw_status = init_branch_status(sw_status_table.eid);
    cost_sw = 1.1;
    cost_load_unsupply=20;
    Objective = cost_sw * sum(abs(sw_status(sw_status_table.varid) - init_sw_status)) + (- 1.0 * cost_load_unsupply * sum(is_restored .* bus_ind.Pd));

    %solve
    sol = optimize([operation_constraints, topology_constraints, powerflow_constraints], Objective, sdpsettings('verbose', 2));

    if sol.problem == 0
        solution.success = true;
        solution.supply = sum(value(is_restored) .* bus_ind.Pd);
        solution.Pij = value(Pij);
        solution.Pi = value(Pi);
        solution.Psub = value(Psub);
        solution.restored = value(is_restored);

        % 还原支路信息
        br_ind.status= true(Nbr,1);
        sol_sw=value(sw_status);
        br_ind.status(sw_status_table.eid) = sol_sw(sw_status_table.varid);
        solution.branch_status = br_ind(:, {'name', 'status'});
        solution.changed_switch = sum(abs(solution.branch_status.status - reshape(init_branch_status, [], 1)));

        % 动作参数
        solution.cost_sw = 1.1;
        solution.cost_load_unsupply = 20;
    else
        solution.success = false;
    end
end